##Reproduces all tables and figures from main text

rm(list=ls())
library(bucky)
library(readstata13)
library(dplyr)
library(tidyr)
library(ggplot2)
library(data.table)
library(xtable)
library(marginaleffects)

completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}

location <- "/Users/dillonlaaker/Dropbox/shocks_replication/"
setwd(paste(location, "data/final", sep = ""))


load("ess_main.Rdata")

####################################
##############FIGURE 1##############
####################################
avg_1825a <- ess %>% filter(age>25) %>%
	group_by(shock4_1825) %>%
	summarise(Mean = weighted.mean(immigration, pspwght, na.rm = TRUE)) %>%
	as.data.frame()
avg_1825a$country <- c("Full Sample", "Full Sample")
avg_1825l <- ess %>% filter(age>25 & datalim==0) %>%
	group_by(shock4_1825) %>%
	summarise(Mean = weighted.mean(immigration, pspwght, na.rm = TRUE)) %>%
	as.data.frame()
avg_1825l$country <- c("Restricted Sample", "Restricted Sample")
	
figure1 <- rbind(avg_1825a, avg_1825l)
figure1$shock4_1825 <- rep(c("No Shock", "Shock"), 2)
colnames(figure1) <- c("Legend", "y", "z")

setwd(paste(location, "paper/figures", sep = ""))
pdf("figure1.pdf", height=3, width=6)
par(mar=c(5,5,0,.25))
dodge <- position_dodge(width = .5)
ggplot(figure1, aes(z, y, fill=Legend)) + 
  ylim(0,2.5) +
  xlab("") + ylab("Anti-Immigration Attitudes")  +
  geom_bar(width = .5, stat = "identity", position=dodge) +
  theme_bw(base_size=10) +
  theme(axis.text = element_text(size=10)) + 
  theme(axis.title.y=element_text(size=10)) +
  scale_fill_manual(values=c("#999999", "#333333")) + 
  theme(panel.grid.minor.x=element_blank(),
  		panel.grid.major.x=element_blank()) + 
  theme(legend.title = element_blank(), legend.position="top") +
  annotate("text", x=0.88, y=round(figure1[1,2], 2) + 0.08, label= as.character(round(figure1[1,2], 2))) + 
  annotate("text", x=1.13, y=round(figure1[2,2], 2) + 0.08, label= as.character(round(figure1[2,2], 2))) + 
  annotate("text", x=1.88, y=round(figure1[3,2], 2) + 0.08, label= as.character(round(figure1[3,2], 2))) + 
  annotate("text", x=2.13, y=round(figure1[4,2], 2) + 0.08, label= as.character(round(figure1[4,2], 2)))
dev.off()

###################################
##############TABLE 2##############
###################################

baseline <- " + secondary + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + aveedu_1017 + eduq_1017 + factor(cohort) + age + I(age^2) + factor(country) + factor(essround)"
##use throughout paper

treatment <- c("shock4_1825", "log(shock4_1825comb+1)")
mod <- paste0("immigration ~ ", treatment, baseline, sep="")
label <- c("Economic Shock", "Economic Shock Count")
table <- data.frame()
counter <- 1

for(d in 1:2){for(i in 1:2){
	model <- lm(mod[i], data=subset(ess, age > 25 & datalim<d), weights=pspwght)
	model <- robustify(model, cluster=c)
	table[counter, 1] <- label[i]
	table[counter, 2] <- summary(model)$coefficients[2,1]
	table[counter, 3] <- summary(model)$coefficients[2,2]
	table[counter, 4] <- nobs(model)
	table[counter, 5] <- ifelse(d==1, "Restricted", "Full")
	print(counter)
	counter <- counter + 1
	}}
table <- table %>% mutate(lci = V2 - 1.96*V3) %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,6:7)] <- format(round(table[, c(2:3,6:7)], 3), nsmall = 3)
table <- table %>% mutate(V3 = paste0("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:5)]
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/table2.tex", sep = ""))


####################################
##############FIGURE 2##############
####################################

y <- c("29", "1017", "1825", "2633", "3441", "4249", "5057")
t1 <- paste0("shock4_", y, sep="")
t2 <- paste0("log(shock4_", y, "comb+1)", sep="")
t <- c(t1, t2)
a <- c(17, 17, 25, 25, 33, 41, 49, 17, 17, 25, 25, 33, 41, 49)

m <- c("aveedu_", "eduq_")
yy <- c("01", "29", "1017", "1825", "2633", "3441", "4249")
v <- NULL
for(i in 1:length(yy)){
	p <- yy[i]
	tt <- paste0(m, rep(p, times=3), sep="", collapse=" + ")
	tt <- paste0(" + ", tt, sep="")

	v <- c(v, tt)
}
v <- rep(v, times=2)
s <- rep("+ shock4_1825", times=7)
s <- c(s, rep("+ log(shock4_1825comb + 1)", times=7))
baselinealt <- " + secondary + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + factor(cohort) + age + I(age^2) + factor(country) + factor(essround)"

counter <- 1
table <- data.frame()

for(d in 1:2){for(i in 1:length(t)){
	mod <- paste0("immigration ~ ", t[i], v[i], baselinealt, sep="")
	model <- lm(mod, data=subset(ess, age > a[i] & datalim<d), weights=pspwght)
	model <- robustify(model, cluster=c)
	table[counter, 1] <- summary(model)$coefficients[2,1]
	table[counter, 2] <- summary(model)$coefficients[2,1] - 1.96*summary(model)$coefficients[2,2]
	table[counter, 3] <- summary(model)$coefficients[2,1] + 1.96*summary(model)$coefficients[2,2]
	table[counter, 4] <- summary(model)$coefficients[2,1] - 1.645*summary(model)$coefficients[2,2]
	table[counter, 5] <- summary(model)$coefficients[2,1] + 1.645*summary(model)$coefficients[2,2]
	table[counter, 6] <- t[i]
	table[counter, 7] <- ifelse(d==1, "Restricted", "Full")
	print(counter)
	counter <- counter + 1
	}}

dr <- table[1:7,]
cr <- table[8:14,]
df <- table[15:21,]
cf <- table[22:28,]ss

setwd(paste(location, "paper/figures", sep = ""))
pdf("figure2.pdf", height=6, width=12)
par(mfrow=c(1,2))
par(mgp=c(3,1,0))
par(mar=c(5,5,0,.25))
counter <- 0
plot(x = NA, xlim = c(-0.12, 0.19), ylim = c(-0.25,6), bty = "n", yaxt = "n", xaxt="n", ylab = "", xlab = "")
axis(side=2, at=c(0,1,2,3,4,5,6), labels=c("2-9","10-17","18-25","26-33","34-41","42-49","50-57"), pos=-0.13, lty=1, las=2)
axis(side=1, at=c(-0.12, -0.06, 0.00, 0.06, 0.12, 0.18), lty=1, las=0)
abline(v = 0, col = "grey", lty = 2)
mtext(side=2, text="Model Specification", line=4, cex=1.25)
mtext(side=1, text="Coefficient on Economic Shock", line=3, cex=1.25)
for(i in 1:nrow(dr)){
	points(x = dr[i, 1], y = counter, pch = 19, cex=1.25)
	segments(x0 = dr[i, 2], x1 = dr[i, 3], y0 = counter, cex=1.25)
	segments(x0 = dr[i, 4], x1 = dr[i, 5], y0 = counter, lwd=2)
	points(x = df[i, 1], y = counter - 0.25, pch = 19, cex=1.25, col="gray")
	segments(x0 = df[i, 2], x1 = df[i, 3], y0 = counter - 0.25, cex=1.25, col="gray")
	segments(x0 = df[i, 4], x1 = df[i, 5], y0 = counter - 0.25, lwd=2, col="gray")
	counter <- counter + 1	
}

par(mgp=c(3,1,0)) 
par(mar=c(5,.25,0,5))
counter <- 0
plot(x = NA, xlim = c(-0.12, 0.19), ylim = c(-.25,6), bty = "n", yaxt = "n", xaxt="n", ylab = "", xlab = "")
axis(side=4, at=c(0,1,2,3,4,5,6), labels=c("2-9","10-17","18-25","26-33","34-41","42-49","50-57"), pos=0.19, lty=1, las=2)
mtext(side=1, text="Coefficient on Economic Shock Count", line=3, cex=1.25)
axis(side=1, at=c(-0.12, -0.06, 0.00, 0.06, 0.12, 0.18), lty=1, las=0)
abline(v = 0, col = "grey", lty = 2)
for(i in 1:nrow(cr)){
	points(x = cr[i, 1], y = counter, pch = 19, cex=1.25)
	segments(x0 = cr[i, 2], x1 = cr[i, 3], y0 = counter, cex=1.25)
	segments(x0 = cr[i, 4], x1 = cr[i, 5], y0 = counter, lwd=2)
	points(x = cf[i, 1], y = counter - 0.25, pch = 19, cex=1.25, col="gray")
	segments(x0 = cf[i, 2], x1 = cf[i, 3], y0 = counter - 0.25, cex=1.25, col="gray")
	segments(x0 = cf[i, 4], x1 = cf[i, 5], y0 = counter - 0.25, lwd=2, col="gray")
	counter <- counter + 1	
}
dev.off()

##################################
#############Figure 3#############
##################################

t <- c(1:6)
t <- paste0("shock", t, "_1825", sep="")
counter <- 1
table <- data.frame()

for(d in 1:2){for(i in 1:length(t)){
	mod <- paste0("immigration ~ ", t[i], baseline, sep="")
	model <- lm(mod, data=subset(ess, age > 25 & datalim<d), weights=pspwght)
	model <- robustify(model, cluster=c)
	table[counter, 1] <- summary(model)$coefficients[2,1]
	table[counter, 2] <- summary(model)$coefficients[2,1] - 1.96*summary(model)$coefficients[2,2]
	table[counter, 3] <- summary(model)$coefficients[2,1] + 1.96*summary(model)$coefficients[2,2]
	table[counter, 4] <- summary(model)$coefficients[2,1] - 1.645*summary(model)$coefficients[2,2]
	table[counter, 5] <- summary(model)$coefficients[2,1] + 1.645*summary(model)$coefficients[2,2]
	table[counter, 6] <- t[i]
	table[counter, 7] <- ifelse(d==1, "Restricted", "Full")
	print(counter)
	counter <- counter + 1
	}}

tab1 <- table[1:6,]
tab2 <- table[7:12,]

setwd(paste(location, "paper/figures", sep = ""))
pdf("figure3.pdf", height=5, width=10)
par(mfrow=c(1,1))
par(mgp=c(3,1,0))
par(mar=c(5,6,0,.25))
plot(x = NA, xlim = c(-0.05, 0.20), ylim = c(-.25,5), bty = "n", yaxt = "n", xaxt="n", ylab = "", xlab = "")
axis(side=2, at=c(0,1,2,3,4,5), labels=c("1 Percent","2 Percent","3 Percent","4 Percent","5 Percent","6 Percent"), pos=-0.055, lty=1, las=1)
axis(side=1, at=c(-0.05, 0.00, 0.05, 0.10, 0.15, 0.20), lty=1, las=1)
abline(v = 0, col = "grey", lty = 2)
mtext(side=2, text="Contraction Threshold", line=5, cex=1.25)
mtext(side=1, text="Coefficient on Economic Shock", line=3, cex=1.25)
counter <- 0
for(i in 1:nrow(table)){
	points(x = tab1[i,1], y = counter+0.10, pch = 19, cex=1.25,)
	segments(x0 = tab1[i,2], x1 = tab1[i,3], y0 = counter+0.10, cex=1.25)
	segments(x0 = tab1[i,4], x1 = tab1[i,5], y0 = counter+0.10, lwd=2)
	points(x = tab2[i,1], y = counter-0.10, pch = 19, cex=1.25, col="gray")
	segments(x0 = tab2[i,2], x1 = tab2[i,3], y0 = counter-0.10, cex=1.25, col="gray")
	segments(x0 = tab2[i,4], x1 = tab2[i,5], y0 = counter-0.10, lwd=2, col="gray")
	counter <- counter + 1
}
dev.off()

###################################
##############TABLE 3##############
###################################

t <- c("shockcont_1825", "unemployment_1825")
counter <- 1
table <- data.frame()

for(i in 1:length(t)){for(d in 1:2){
	mod <- paste0("immigration ~ ", t[i], baseline, sep="")
	model <- lm(mod, data=subset(ess, age > 25 & datalim<d), weights=pspwght)
	model <- robustify(model, cluster=c)
	table[counter, 1] <- ifelse(t[i]=="shockcont_1825", "Economic Shock Cont.", "Unemployment Rate")
	table[counter, 2] <- summary(model)$coefficients[2,1]
	table[counter, 3] <- summary(model)$coefficients[2,2]
	table[counter, 4] <- nobs(model)
	table[counter, 5] <- ifelse(d==1, "Restricted", "Full")
	print(counter)
	counter <- counter + 1
	}}

table <- table %>% mutate(lci = V2 - 1.96*V3) %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,6:7)] <- format(round(table[, c(2:3,6:7)], 3), nsmall = 3)
table$lci <- gsub(" ", "", table$lci)
table <- table %>% mutate(V3 = paste0("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:5)]
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/table3.tex", sep = ""))


##################################
#############Figure 4#############
##################################

int <- c(rep("secondary", times=4), rep("skill", times=6), rep("pedu", times=10))
sample <- c(rep("restricted", times=2), rep("full", times=2), rep("restricted", times=3), rep("full", times=3), rep("restricted", times=5), rep("full", times=5))
value <- c("12 years or more", "Less than 12 years", "12 years or more", "Less than 12 years", "High skilled", "Mod. skilled", "Low skilled", "High skilled", "Mod. skilled", "Low skilled", "Less educated", "2", "4", "6", "More educated", "Less educated", "2", "4", "6", "More educated")
est <- rep(NA, times=20)
se <- rep(NA, times=20)
table <- data.frame(int, sample, value, est, se)

##PlotA
model <- lm(immigration ~ shock4_1825*secondary + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + aveedu_1017 + eduq_1017 + factor(cohort) + age + I(age^2) + factor(country) + factor(essround), data=subset(ess, age > 25 & datalim<1), weights=pspwght)
model <- robustify(model, cluster=c)
model <- slopes(model, variables="shock4_1825", by="secondary")
table[1,4] <- model[1,4]
table[2,4] <- model[2,4]
table[1,5] <- model[1,5]
table[2,5] <- model[2,5]

model <- lm(immigration ~ shock4_1825*secondary + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + aveedu_1017 + eduq_1017 + factor(cohort) + age + I(age^2) + factor(country) + factor(essround), data=subset(ess, age > 25 & datalim<2), weights=pspwght)
model <- robustify(model, cluster=c)
model <- slopes(model, variables="shock4_1825", by="secondary")
table[3,4] <- model[1,4]
table[4,4] <- model[2,4]
table[3,5] <- model[1,5]
table[4,5] <- model[2,5]

##PlotB
model <- lm(immigration ~ shock4_1825*factor(foccupation14) + secondary + factor(moccupation14) + female + minority + pcitizen + peduc + aveedu_1017 + eduq_1017 + factor(cohort) + age + I(age^2) + factor(country) + factor(essround), data=subset(ess, age > 25 & datalim<1), weights=pspwght)
model <- robustify(model, cluster=c)
model <- slopes(model, variables="shock4_1825", by="foccupation14")
table[5,4] <- model[1,4]
table[6,4] <- model[3,4]
table[7,4] <- model[2,4]
table[5,5] <- model[1,5]
table[6,5] <- model[3,5]
table[7,5] <- model[2,5]

model <- lm(immigration ~ shock4_1825*factor(foccupation14) + secondary + factor(moccupation14) + female + minority + pcitizen + peduc + aveedu_1017 + eduq_1017 + factor(cohort) + age + I(age^2) + factor(country) + factor(essround), data=subset(ess, age > 25 & datalim<2), weights=pspwght)
model <- robustify(model, cluster=c)
model <- slopes(model, variables="shock4_1825", by="foccupation14")
table[8,4] <- model[1,4]
table[9,4] <- model[3,4]
table[10,4] <- model[2,4]
table[8,5] <- model[1,5]
table[9,5] <- model[3,5]
table[10,5] <- model[2,5]

##PlotC
model <- lm(immigration ~ shock4_1825*peduc + secondary + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + aveedu_1017 + eduq_1017 + factor(cohort) + age + I(age^2) + factor(country) + factor(essround), data=subset(ess, age > 25 & datalim<1), weights=pspwght)
model <- robustify(model, cluster=c)
model <- slopes(model, variables="shock4_1825", by="peduc")
table[11,4] <- model[7,4] ## peduc==0
table[12,4] <- model[6,4] ## peduc==2
table[13,4] <- model[3,4] ## peduc==4
table[14,4] <- model[1,4] ## peduc==6
table[15,4] <- model[8,4] ## peduc==8
table[11,5] <- model[7,5] ## peduc==0
table[12,5] <- model[6,5] ## peduc==2
table[13,5] <- model[3,5] ## peduc==4
table[14,5] <- model[1,5] ## peduc==6
table[15,5] <- model[8,5] ## peduc==8

model <- lm(immigration ~ shock4_1825*peduc + secondary + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + aveedu_1017 + eduq_1017 + factor(cohort) + age + I(age^2) + factor(country) + factor(essround), data=subset(ess, age > 25 & datalim<2), weights=pspwght)
model <- robustify(model, cluster=c)
model <- slopes(model, variables="shock4_1825", by="peduc")
table[16,4] <- model[7,4] ## peduc==0
table[17,4] <- model[6,4] ## peduc==2
table[18,4] <- model[3,4] ## peduc==4
table[19,4] <- model[1,4] ## peduc==6
table[20,4] <- model[8,4] ## peduc==8
table[16,5] <- model[7,5] ## peduc==0
table[17,5] <- model[6,5] ## peduc==2
table[18,5] <- model[3,5] ## peduc==4
table[19,5] <- model[1,5] ## peduc==6
table[20,5] <- model[8,5] ## peduc==8

table$lci95 <- table[,4] - table[,5] * 1.96
table$uci95 <- table[,4] + table[,5] * 1.96
table$lci90 <- table[,4] - table[,5] * 1.645
table$uci90 <- table[,4] + table[,5] * 1.645

hs <- table[c(1:4),]
occup <- table[c(5:10),]
peduc1 <- table[c(11:15),]
peduc2 <- table[c(16:20),]

setwd(paste(location, "paper/figures", sep = ""))
pdf("figure4.pdf", height=5, width=10)
par(mfrow=c(1,3))
par(mar=c(5,5,0,.25))
plot(x = NA, ylim = c(-0.10, 0.27), xlim = c(-.05,1.05), bty = "n", yaxt = "n", xaxt="n", ylab = "", xlab = "")
axis(side=1, at=c(0.10,0.90), labels=c("Less than\n12 years", "12 years\nor more"), pos=-0.105, lty=1, las=1, cex.axis=1.15, mgp = c(3,2,0))
axis(side=2, at=c(-0.10, -0.05, 0.00, 0.05, 0.10, 0.15, 0.20, 0.25), lty=1, las=0, cex.axis=1.15)
abline(h = 0, col = "grey", lty = 2)
mtext(side=1, text="Respondent's Education", line=3.5, cex=1.25)
mtext(side=2, text="Marginal Effect of Economic Shock", line=3, cex=1.25)
counter <- 0.10
for(i in 2:1){
  points(y = hs[i,4], x = counter+0.025, pch = 19, cex=1.25,)
  segments(y0 = hs[i,6], y1 = hs[i,7], x0 = counter+0.025, cex=1.25)
  segments(y0 = hs[i,8], y1 = hs[i,9], x0 = counter+0.025, lwd=2)
  points(y = hs[i+2,4], x = counter-0.025, pch = 19, cex=1.25, col="gray")
  segments(y0 = hs[i+2,6], y1 = hs[i+2,7], x0 = counter-0.025, cex=1.25, col="gray")
  segments(y0 = hs[i+2,8], y1 = hs[i+2,9], x0 = counter-0.025, lwd=2, col="gray")
  counter <- counter + 0.80
}

par(mar=c(5,4,0,.25))
plot(x = NA, ylim = c(-0.10, 0.27), xlim = c(-.05,2.05), bty = "n", yaxt = "n", xaxt="n", ylab = "", xlab = "")
axis(side=1, at=c(0.20,1,1.80), labels=c("High\nskilled", "Mod.\nskilled", "Low\nskilled"), pos=-0.105, lty=1, las=1, cex.axis=1.15, mgp = c(3,2,0))
axis(side=2, at=c(-0.10, -0.05, 0.00, 0.05, 0.10, 0.15, 0.20, 0.25), lty=1, las=0, cex.axis=1.15)
abline(h = 0, col = "grey", lty = 2)
mtext(side=1, text="Parental occupation", line=3.5, cex=1.25)
counter <- 0.20
for(i in 1:3){
  points(y = occup[i,4], x = counter+0.05, pch = 19, cex=1.25,)
  segments(y0 = occup[i,6], y1 = occup[i,7], x0 = counter+0.05, cex=1.25)
  segments(y0 = occup[i,8], y1 = occup[i,9], x0 = counter+0.05, lwd=2)
  points(y = occup[i+3,4], x = counter-0.05, pch = 19, cex=1.25, col="gray")
  segments(y0 = occup[i+3,6], y1 = occup[i+3,7], x0 = counter-0.05, cex=1.25, col="gray")
  segments(y0 = occup[i+3,8], y1 = occup[i+3,9], x0 = counter-0.05, lwd=2, col="gray")
  counter <- counter + 0.80
}

par(mar=c(5,4,0,.25))
plot(x = NA, ylim = c(-0.10, 0.27), xlim = c(-1,9), bty = "n", yaxt = "n", xaxt="n", ylab = "", xlab = "")
axis(side=1, at=c(0,2,4,6,8), labels=c("Less\neducated", "2", "4", "6", "More\neducated"), pos=-0.105, lty=1, las=1, cex.axis=1.15, mgp = c(3,2,0))
axis(side=2, at=c(-0.10, -0.05, 0.00, 0.05, 0.10, 0.15, 0.20, 0.25), lty=1, las=0, cex.axis=1.15)
abline(h = 0, col = "grey", lty = 2)
mtext(side=1, text="Parental education", line=3.5, cex=1.25)
counter <- 0
for(i in 1:nrow(peduc1)){
  points(y = peduc1[i,4], x = counter+0.25, pch = 19, cex=1.25,)
  segments(y0 = peduc1[i,6], y1 = peduc1[i,7], x0 = counter+0.25, cex=1.25)
  segments(y0 = peduc1[i,8], y1 = peduc1[i,9], x0 = counter+0.25, lwd=2)
  points(y = peduc2[i,4], x = counter-0.25, pch = 19, cex=1.25, col="gray")
  segments(y0 = peduc2[i,6], y1 = peduc2[i,7], x0 = counter-0.25, cex=1.25, col="gray")
  segments(y0 = peduc2[i,8], y1 = peduc2[i,9], x0 = counter-0.25, lwd=2, col="gray")
  counter <- counter + 2
}
dev.off()

##################################
#############Figure 5#############
##################################

#takes a long time to run.

ess <- select(ess, immigration, shock4_1825, shock4_1825comb, female, minority, pcitizen, peduc, age, cohort, country, educationy, locationtype, union, foccupation14, moccupation14, lrideology, income, essround, pspwght, c, datalim)
ess$age[ess$age<26] <- NA
ess$datalim[ess$datalim==1] <- NA
ess <- na.omit(ess)

set.seed(12345)
boots <- 1000
acde.boots_edu <- acde.boots_inc <- acde.boots_ideo <- pt.boots <- rep(NA, times = boots)
munis <- unique(ess$c)
lookup <- split(1:nrow(ess), ess$c)
codes <- names(lookup)
for (b in 1:boots) {
	draw <- sample(codes, length(codes), replace = TRUE)
	star <- unlist(lookup[draw])
	ess.star <- ess[star,]
  	Sys.sleep(0.1)
    print(b)
    # update GUI console
    flush.console()
    
    boot.first <- lm(immigration ~ shock4_1825 + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + age + I(age^2) + factor(cohort) + factor(country) + educationy + locationtype + union + lrideology + income + factor(essround), data=ess.star, weights=pspwght)
    pt.boots[b] <- coef(boot.first)["shock4_1825"]
    
    boot.direct <- lm(I(immigration - coef(boot.first)["educationy"]*educationy) ~ shock4_1825 + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + age + I(age^2) + factor(cohort) + factor(essround) + factor(country), data=ess.star, weights=pspwght)
    acde.boots_edu[b] <- coef(boot.direct)["shock4_1825"]
    
    boot.direct <- lm(I(immigration - coef(boot.first)["income"]*income) ~ shock4_1825 + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + age + I(age^2) + factor(cohort) + factor(essround) + factor(country), data=ess.star, weights=pspwght) 
    acde.boots_inc[b] <- coef(boot.direct)["shock4_1825"]

    boot.direct <- lm(I(immigration - coef(boot.first)["lrideology"]*lrideology) ~ shock4_1825 + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + age + I(age^2) + factor(cohort) + factor(essround) + factor(country), data=ess.star, weights=pspwght)
	  acde.boots_ideo[b] <- coef(boot.direct)["shock4_1825"]
}

setwd(paste(location, "analysis/acde", sep = ""))
write.csv(acde.boots_edu, file="acde_education_restricted.csv", row.names=F)
write.csv(acde.boots_inc, file="acde_income_restricted.csv", row.names=F)
write.csv(acde.boots_ideo, file="acde_ideology_restricted.csv", row.names=F)
write.csv(pt.boots, file="ptboots_restricted.csv", row.names=F)


setwd(paste(location, "data/final", sep = ""))
load("ess_main.Rdata")
ess <- select(ess, immigration, shock4_1825, shock4_1825comb, female, minority, pcitizen, peduc, age, cohort, country, educationy, locationtype, union, foccupation14, moccupation14, lrideology, rti, offshorable, manufacture, income, essround, pspwght, c, datalim)
ess$age[ess$age<26] <- NA
ess <- na.omit(ess)

set.seed(12345)
boots <- 1000
acde.boots_edu <- acde.boots_inc <- acde.boots_ideo <- pt.boots <- rep(NA, times = boots)
munis <- unique(ess$c)
lookup <- split(1:nrow(ess), ess$c)
codes <- names(lookup)
for (b in 1:boots) {
  draw <- sample(codes, length(codes), replace = TRUE)
  star <- unlist(lookup[draw])
  ess.star <- ess[star,]
  ##ess.star <- ess[sample(1:nrow(ess), replace = TRUE),]
  Sys.sleep(0.1)
  print(b)
  # update GUI console
  flush.console()
  
  boot.first <- lm(immigration ~ shock4_1825 + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + age + I(age^2) + factor(cohort) + factor(country) + educationy + locationtype + union + lrideology + income + factor(essround), data=ess.star, weights=pspwght)
  pt.boots[b] <- coef(boot.first)["shock4_1825"]
  
  boot.direct <- lm(I(immigration - coef(boot.first)["educationy"]*educationy) ~ shock4_1825 + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + age + I(age^2) + factor(cohort) + factor(essround) + factor(country), data=ess.star, weights=pspwght)
  acde.boots_edu[b] <- coef(boot.direct)["shock4_1825"]
  
  boot.direct <- lm(I(immigration - coef(boot.first)["income"]*income) ~ shock4_1825 + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + age + I(age^2) + factor(cohort) + factor(essround) + factor(country), data=ess.star, weights=pspwght) 
  acde.boots_inc[b] <- coef(boot.direct)["shock4_1825"]
  
  boot.direct <- lm(I(immigration - coef(boot.first)["lrideology"]*lrideology) ~ shock4_1825 + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + age + I(age^2) + factor(cohort) + factor(essround) + factor(country), data=ess.star, weights=pspwght)
  acde.boots_ideo[b] <- coef(boot.direct)["shock4_1825"]
}

setwd(paste(location, "code/analysis/acde", sep = ""))
write.csv(acde.boots_edu, file="acde_education_full.csv", row.names=F)
write.csv(acde.boots_inc, file="acde_income_full.csv", row.names=F)
write.csv(acde.boots_ideo, file="acde_ideology_full.csv", row.names=F)
write.csv(pt.boots, file="ptboots_full.csv", row.names=F)

files <- c("acde_income_full.csv", "acde_ideology_full.csv", "acde_education_full.csv", "acde_income_restricted.csv", "acde_ideology_restricted.csv", "acde_education_restricted.csv")
table <- data.frame()
n <- c("Net Income (Full)", "Net Ideology (Full)", "Net Education (Full)", "Net Income (Restricted)", "Net Ideology (Restricted)", "Net Education (Restricted)")

for(i in 1:length(files)){
	data <- read.csv(files[i])
	table[i,1] <- mean(data$x)
	table[i,2] <- quantile(data, probs = c(0.025), na.rm=TRUE)
	table[i,3] <- quantile(data, probs = c(0.975), na.rm=TRUE)
	table[i,4] <- n[i]	
}

table1 <- table[c(1:3),]
table1 <- rbind(c(0.082,0.040,0.123,"Original (Full)"), table1)
table2 <- table[c(4:6),]
table2 <- rbind(c(0.118, 0.073, 0.164,"Original (Restricted)"), table2)


setwd(paste(location, "paper/figures", sep = ""))
pdf("figure5.pdf", height=4, width=9) 
par(mfrow=c(1,1))
par(mgp=c(3,1,0))
par(mar=c(5,1,1,0))
plot(x = NA, xlim = c(-0.05, 0.205), ylim = c(-.5,11.5), bty = "n", yaxt = "n", xaxt="n", ylab = "", xlab = "")
mtext("ACDE of Economic Shock", side=1, line=2.5, cex=1.05)
axis(side=1, at=c(-0.05, 0.00, 0.05, 0.10, 0.15, 0.20), lty=1, las=0)
abline(v = 0, col = "grey", lty = 2)
counter <- 11
for(i in 1:length(table1)){
	points(x = table2[i,1], y = counter, pch = 19, cex=0.80)
	segments(x0 = as.numeric(table2[i,2]), x1 = as.numeric(table2[i,3]), y0 = counter, cex=0.80)
	text(x = table2[i,1], y = counter+.5, table2[i,4], cex=0.80)
	points(x = table1[i,1], y = counter-6, pch = 19, cex=0.80)
	segments(x0 = as.numeric(table1[i,2]), x1 = as.numeric(table1[i,3]), y0 = counter-6, cex=0.80)
	text(x = table1[i,1], y = counter-6+.5, table1[i,4], cex=0.80)
	counter <- counter - 1.5
	
}
dev.off()





